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We present the results of molecular dynamic simulations 
of a two-dimensional vortex array driven by a uniform cur- 
rent through random pinning centers at zero temperature. 
We identify two types of flow of the driven array near the 
depinning threshold. For weak disorder the flux array con- 
tains few dislocation and moves via correlated displacements 
of patches of vortices in a crinkle motion. As the disorder 
strength increases, we observe a crossover to a spatially inho- 
mogeneous regime of plastic flow, with a very defective vortex 
array and a channel-like structure of the flowing regions. The 
two regimes are characterized by qualitatively different spa- 
tial distribution of vortex velocities. In the crinkle regime the 
distribution of vortex velocities near threshold has a single 
maximum that shifts to larger velocities as the driving force 
is increased. In the plastic regime the distribution of vortex 
velocities near threshold has a clear bimodal structure that 
persists upon time-averaging the individual velocities. The 
bimodal structure of the velocity distribution reflects the co- 
existence of pinned and flowing regions and is proposed as a 
quantitative signature of plastic flow. 
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I. INTRODUCTION 

The problem of nonlinear collective transport through 
random media has attracted much theoretical and exper- 
imental attention due to the interesting spatio-temporal 
phenomena that arise in a variety of physical systems 
from the competition between interactions and disorder. 
In particular, the dynamics of driven elastic media that 
are distorted by disorder, but cannot "break" , has been 
studied extensively, both theoretically and numerically. 
At zero temperature these systems exhibit a sharp de- 
pinning transition from a pinned state below a critical 
driving force Fc to a sliding state above Fc- The transi- 
tion can be described as a critical phenomenon in terms of 
scaling laws and critical exponents. The elastic medium 
model can be used to describe the dynamics of weakly 
pinned Abrikosov flux lattices , fronts of wetting flu- 
ids invading porous media |^ and charge density waves 
(CDW's) in anisotropic conductors over a wide range 
of length scales. It is, however, expected to eventu- 
ally break down (particularly in lower dimensionality) 
at short length scales since it yields unphysical regions 
of unbounded strains 1^ . In addition, the elastic model 
is inadequate for many physical systems with strong dis- 
order that exhibit a spatially inhomogeneous plastic re- 



sponse without long-wavelength elastic restoring forces. 
These include strongly pinned flux lattices |^-^ , invasion 
of nonwetting fluid in porous media |Q , Wigner solids in 
2DEG and fluid flow down a rough incline Q. In these 
systems the competition between drive and disorder gen- 
erates defects (dislocations, phase slips) in the medium 
that can qualitatively change the dynamics pO|-p^ . Col- 
lective transport in the presence of topological defects is 
still poorly understood. 

Magnetic flux arrays in type-II superconductors are an 
ideal system for investigating nonlinear collective trans- 
port since by changing the applied magnetic field one 
can tune the strength of the intervortex interaction and 
observe a crossover from a regime of weak pinning, well 
described by collective flux creep theories, to a regime 
of strong pinning with spatially inhomogeneous flow [0 . 
Evidence for this comes from early simulations of two- 
dimensional flux lattices by Jensen et al. |^ and by Shi 
and Berlinsky In addition, a variety of transport phe- 
nomena observed recently in superconductors has been 
attributed to the inhomogeneous plastic response of the 
flux array, including a nonmonotonic dependence of the 
critical current on temperature just below the flux lattice 
melting point (peak effect), an unusual current and field 
dependence of 1// broadband noise and fingerprint phe- 
nomena Isyi^^jl^. On the other hand, the experiments 
only probe flux motion indirectly through transport mea- 
surements. For this reason their interpretation is difficult 
and still controversial. Numerical work can therefore be 
very valuable to gain insight into this complex problem 
and to serve as a guide for future theoretical work. 

In this paper we report on simulations of the dynam- 
ics of a two-dimensional Abrikosov fiux array driven by a 
uniform current through random pinning centers at zero 
temperature. The focus of our work, which distinguishes 
it from previous numerical work on the same 0,^,^,^ 
or closely related [qI JTtIJi^ systems, is on identifying var- 
ious types of flow and establishing a connection between 
the type of flow or response ("elastic" versus "plastic") 
and the presence of flux lattice defects and the shape of 
the macroscopic response as embodied, for instance, in 
the V-I characteristics. This will be very useful for the 
interpretation of experiments. While most of the results 
presented here are somewhat qualitative, our long-term 
objective is to carry out simulations for realistic param- 
eter values that will allow a detailed comparison with 
experiments. It is well known that short-wavelength de- 
fects, such as dislocations, play a more important role in 
two, rather than in three, dimensions pM- For this rea- 
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son many of our results will not apply directly to three- 
dimensional flux arrays. On the other hand, the study 
of two-dimensional systems is valuable both because of 
intrinsic interest and because in many experimental sit- 
uations the flux arrays can effectively be modeled as two 
dimensional. Thermal fluctuations are generally impor- 
tant in flux flow experiments and purely dynamical phe- 
nomena associated with the current-induced depinning of 
the vortices cannot be dissociated from thermal-induced 
softening of the lattice. In this paper we speciflcally con- 
sider the flux array dynamics at T = with the objective 
of disentangling these two types of effects. 

Below we present the results of simulations of a driven 
flux array for both a low density and a very high den- 
sity of point pinning centers. In both cases we identify 
two types of response or flow of the flux array near the 
depinning threshold and a crossover from one type to 
another as the disorder strength is increased. For suf- 
ficiently weak disorder the flux array contains very few 
defects and moves via correlated displacements of patches 
of vortices. The dynamics is similar to that observed by 
Hu and Westervelt |Q in magnetic bubble arrays. Fol- 
lowing these authors' suggestion, we refer to this type 
of response as crinkle motion. For stronger disorder the 
response near threshold is plastic, with vortices flowing 
around pinned regions. The flux lattice is very defective 
and we observe channels of liquid-like vortex array flow- 
ing around solid-like pinned regions. The topology of the 
channels is not, however, fixed in time. Channels open 
and close continuously as the flux array is driven over the 
impurities and all vortices participate in the motion at 
one time or another near threshold. As the disorder is 
further increased the individual channels become longer 
lived and for very strong disorder we observe a filamen- 
tary structure with a fraction of vortices that never move 
on the time scale of the simulation. To characterize the 
different regimes we have studied in detail the spatial 
distribution of vortex velocities. In the crinkle regime 
the distribution of vortex velocity near threshold shows 
a single maximum corresponding roughly to the average 
velocity of the array, though at any time, some vortices 
are moving with velocity significantly greater than the 
average value. The plastic fiow regime is characterized 
by bimodal velocity distributions near threshold, indi- 
cating that the velocity is spatially inhomogeneous, with 
both pinned and fiowing regions. We discuss the corre- 
lation between these qualitative features of the velocity 
distribution and the shape of the macroscopic V-I char- 
acteristic and suggest that the shape of the velocity dis- 
tribution may be used for a crude classification of the 
type of response. 



II. THE MODEL 

The specific model considered here is essentially the 
same as that studied in earlier simulations by Jensen et 
al. 01, by Shi and Berlinsky |^ and, more recently, by 
Koshelev and Vinokur |TT| . The two-dimensional pan- 
cake vortices are modeled as point-particles with finite- 
range interaction and overdamped dynamics, driven 
through randomly placed pinning centers by a uniform 
force F proportional to the external current. The equa- 
tions of motion for the two-dimensional vortex positions 
Vi are given by 

ijij k=i 

(1) 

where {R/c} denote the random positions of the Np pin- 
ning centers and Ny is the total number of vortices. Here, 
7i is the friction coefficient of a single vortex, which will 
be incorporated in our unit of time. The repulsive inter- 
vortex interaction has finite range Rc and yields a force 
— VV^(r) = /^e"''/^" (l — r/i?c)r on the zth vortex, with 
Rv < Rc- The results presented below have been ob- 
tained with — Rc- The second term on the right 
hand side of Eq. (|l]) is the attractive pinning force of 
range i?p, given by -Wp{r) = - r'^/RD^r/Rp. 

The range Ry of the intervortex repulsion is of the order 
of the superconductor penetration length. A, while the 
range Rp of the pinning potential is of the order of the 
superconductor coherence length, ^. In the absence of 
pinning the flux array forms a stable triangular lattice of 
lattice constant ag. 

One of the difficulties in carrying out a detailed nu- 
merical study of the dynamical response of this model 
system is the large number of parameters to be varied. 
In the following we have chosen the range Rp of the pin- 
ning potential as our unit of length and the maximum 
intervortex force as our unit of force. In all cases the 
vortex lattice is rather dense, with UyR^ ~ 8 — 9, where 
TLy = l/a§ is the areal density of vortices, and soft, with 
C66 ~ 0.271 — 0.334, where cge is the shear modulus of the 
clean vortex lattice, in the absence of disorder. The pin- 
ning centers are modeled as point pins, in the sense that 
ao,Rv >> Rp. We have considered sets of parameters 
corresponding to two rather different density of pins: (i) 
a dense array of overlapping pins, with Np/N^ ~ 133 and 
UpRp « 2.8 with Tip the areal density of pins, and (ii) a di- 
lute array of nonoverlapping pins, with Np/Ny = 0.5 and 
UpRp ~ 0.046. The specific parameter values used are 
given below. In both cases we have varied the strength 
of the disorder by varying the maximum pinning force, 
/p. 

The mean motion of the flux array is described by the 
drift velocity in the direction F of the driving force, given 
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by 



1 



^ i—1 



(2) 



The angular brackets denote the average over disorder. 
In the numerical calculation we average over impurity 
realizations by performing a time average since as time 
evolves the flux array samples different impurity configu- 
rations. The mean velocity is proportional to the voltage 
V from flux motion, while the driving force F is pro- 
portional to the driving current /. Curves of Vd versus 
F correspond therefore to the V-I characteristics of the 
material. 



III. ELASTIC RESPONSE 

Even in the absence of driving force the random pin- 
ning produces both elastic and plastic (dislocations) de- 
formations of the lattice. If topological defects are ex- 
plicitly forbidden in the model, the distortion induced 
by disorder can be described within elasticity theory. 

Treating the disorder as a perturbation, it has been 
shown that for d < 4 order persists only in region of lin- 
ear size Rc that are pinned collectively [|l] . The pinning 
length Rc can be estimated by an Imry-Ma argument 
p2| by assuming that in the presence of disorder the 
flux array deforms elastically to take advantage of the 
pinning wells. The elastic energy cost associated with 
displacing a region of linear size i? by a distance Rp is 
6Eei{R) ^ cqq{Rp/R)^R'^ in d dimensions, where cge is 
the shear modulus of the flux lattice. The correspond- 
ing pinning energy gain is 5Ep{R) ~ ^/nyT{R/ Rp)'^'^/'^'' 
where T is the variance of the disorder potential arising 
from uncorrelated point pins, with T « np{fpRp)^. For 
dimension d > 4 the elastic energy always exceeds the 
pinning energy at large distances and the ordered state 
is stable for small disorder. For c? < 4 disorder dominates 
beyond the length Rc where the elastic strains induced 
by disorder are of order one, or 5Eci{Rc) 5Ep{Rc). The 
elastic medium is broken up in domains of size Rc, given 
by 



Rc 



RpCee 



(3) 



in d — 2. Alternatively, the collective pinning length 
can be defined following Larkin and Ovchinnikov pl| ] 
as the length scale where the mean square displacement 
< >i/2 induced by disorder is of the order of the 

range of the pinning potential, i.e., < [u{Rc)]'^ >^ Rp. 
The estimate described above neglects logarithmic cor- 
rections and assumes point pins of range Rp << Qq. 

The Larkin-Ovchinnikov collective pinning theory ap- 
plies provided Rc >> a^. In addition, if topological de- 
fects can occur in the lattice, the mean distance between 



such defects must exceed Rc- While our simulations have 
been carried out for parameter values where the above 
inequalities generally do not apply, it is useful to briefly 
summarize some of the properties of driven elastic media 
for comparison. The force needed to depin a region of 
linear size R can be found by equating the energy gain 
due to the external force ~ FRpR^ to the pinning energy 
dEp(R) of the region and is given by 
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RIR 



(4) 



In the weak pinning regime where Larkin domains of size 
Rc are pinned collectively by disorder, the threshold force 
needed to depin the medium can be estimated as the 
force needed to depin a Larkin domain, F™" = F{Rc) ~ 

ni,r/(c66i?p). 

When Rc ^ or ^/V /Rp ^ ceeRp, the collective pin- 
ning theory breaks down and vortices are pinned individ- 
ually. In this strong pinning regime the threshold force 
can be estimated as p^^°^3 _ pl^R ^ ag) ^ ny\/T/R^. 
The disorder-induced displacement of the lattice exceeds 
the range of the pinning potential and the disorder can no 
longer be treated as a perturbation. When this displace- 
ment becomes of order ao, the Fourier components of the 
pinning potential with the periodicity of the underlying 
lattice become dominant and change qualitatively the na- 
ture of the pinning |23 23]. In this case one needs to go 



beyond the simple dimensional estimates just described, 
as discussed recently by Giamarchi and Le Doussal [ p^ . 

The question of when and how, as a function of dis- 
order strength, topological defects proliferate has been 
addressed recently by Gingras and Huse for a ferro- 
magnetic random field XY model. Dislocations allow a 
region of linear size R to better adjust to disorder, yield- 
ing a gain in pinning energy. Gingras and Huse argue 
that a bound for the length scale Rd where dislocations 
proliferate can be obtained by equating the elastic energy 
cost of a dislocation ^ CQ^a^hiR to this pinning energy 
gain. If the pinning energy gain is estimated again as 
5Ep{R) - y/^R/Rp, we obtain Rd ~ (ag/i?^)i?^. No- 
tice, however, that since the disorder-induced displace- 
ments of the lattice in the presence of dislocations exceed 
Rp, the pinning energy no longer grows linearly with dis- 
placement and this estimate is at best a lower bound of 
the actual pinning energy. For this reason all we can re- 
ally infer from this argument is that in the context of 
weak collective pinning Rd > Rc- The focus of our paper 
is not on determining the length Rd, but rather on the 
dynamics of the driven system and on the proliferation or 
healing of dislocations as a result of the competition be- 
tween disorder, drive and interactions. We can estimate 
the force needed to depin and heal dislocations separated 
by a length Ld as Fd ~ F{Ld), where F{L) is given in 
Eq. ^, with the result Fd — UyT / (ceealRp) . We remark 
that this dimensional estimate for Fd is identical to the 
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"crystallization" force Ft of Koshelev and Vinokur. 

We can then distinguish three regions as a function 
of the disorder strength F, as shown in the schematic 
"phase" diagram of Fig. 1. These regions may or 
may not be separated by actual phase transitions. For 
F < (cge^^p)^, Ld > Lc > ao and the pinning is collec- 
tive. The driven medium responds elastically and the rel- 
evant threshold force for depinning is the force ^ F 
needed to depin a Larkin domain. For [cqqR^)^ < F < 
{cQQO^Rpf' , the pinning is strong since Ld > ao > Lc and 
vortices are pinned individually. The threshold force is 
estimated as the force f^*^°"3 ^ needed to depin a 
single vortex. In both these regions the force Fd needed 
to depin and heal dislocations present in the lattice is 
smaller than the threshold force. When F > (ceeOo-^p)^, 
then Ld < ao and the force ~ F exceeds the thresh- 
old force for depinning. In this region one may have a 
scenario of the type proposed by Koshelev and Vinokur 
PI , with a pinned disordered solid for F < F^*'""^, a 
region of plastic flow for p^^°^3 < F < Fd and finally a 
flowing solid with no topological defects for F > Fd- 
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FIG. 1. A schematic "phase diagram" illustrating the var- 
ious regimes in the {F — F) plane. The lines separating the 
various regions represent the estimates of threshold force dis- 
cussed in the text. The boundary between pinned and slid- 
ing regimes is the collective threshold force F^°" ~ F for 
F < Fi ~ (ceeRp)^ and the strong pinning threshold force 
^strong ^ j.^^ F > Fi. For F > r2 ~ {cesaoRpf there is a 
region of plastic flow above the pinned region, separated from 
the sliding solid by the force ~ F. 

The sliding state of an elastic medium driven through a 
random potential can be studied analytically via a high- 
velocity perturbation theory. The perturbation theory 



was introduced by Schmid and Hanger [ p6[ and by Larkin 
in the context of flux lattices and further developed by 
Sneddon, Cross, and Fisher for sliding CDW's. More 
recently Zhu, Littlewood and Millis discussed in detail 
the high- velocity perturbation theory for sliding Wigner 
crystals 

The starting point of the perturbation theory is a 
description of the flux array as an overdamped elastic 
continuum driven by an external force and distorted by 
short-range uncorrelated disorder. The equation of mo- 
tion for the two-dimensional displacement field u(r, t) is 
given by 

jdtu{r, t) = (cii - C66)V(V • u) + ceeV^u -t- Fp(r, t) + n„F, 

(5) 

where 7 is a friction per unit area, related to the single- 
vortex friction coefficient of Eq. |l| by 7 = riy'fi, cn and 
are the comprcssional and shear clastic moduli of the 
two-dimensional flux lattice, and Fp is the pinning force 
per unit area, 

Fp(r,t)--po(r)Vy(r-|-u(r,t)), (6) 

where po{r) — J2n=i '^('^"^^n) the spatially inhomoge- 
neous density of the undistorted lattice, with R° the sites 
of the triangular Abrikosov lattice. The coarse-grained 
quenched pinning potential, V{r) has zero mean and 
short ranged correlations, < V{r)V{r') >= F(r)/(r — r'), 
with /(r) a function that drops rapidly to zero for r >> 
Rp. For simplicity of notation we have neglected in Eq. 
^ the nonlocality of the elastic constants. This can, how- 
ever, be easily incorporated in the perturbation theory. 
The drift velocity is defined here as Vd{F) =< dtU-F >, 
where the brackets denote a spatial average, as well as 
a disorder average. In the absence of disorder Vd{F) = 
F/-f. Treating the disorder as a perturbation relative to 
the external driving force, one can then evaluate correc- 
tions to this uniformly sliding state. The details of the 
calculation are not given here as this follows closely the 
perturbation theory for the Wigner crystal described re- 
cently by Zhu et al. |2^]. Rather than expanding about 
the solution Vd(F) — Ff-f in the absence of disorder, one 
actually constructs a self consistent perturbation theory 
by writing u{r,t) = Vdt + s{r,t), solving for s(r,t) in 
perturbation theory and then requiring < (?ts(r, t) >= 0. 
The lowest non-vanishing correction Svd = f d — F/j to 
the mean velocity is given by, 

'5^'d(F)«-^^F(G)G2(G.F) 

f _dk^ 7vG 

where G are the reciprocal lattice vectors of the trian- 
gular flux lattice and the k-integral is over the first Bril- 
louin zone of size ksz = y/4wri^. The disorder correlator 
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is given by r(G) w T « npifpR^)^, for G < 1/Rp, and 
vanishes rapidly for G >> 1/Rp. It therefore cuts off the 







reciprocal lattice vectors sum at Gmax ~ ^/Rp- For the 
case of point pins (oq >> -Rp), the right hand side of Eq. 

can be evaluated exactly by transforming the wavevec- 
tor sum to an integral. It is, however, more instructive to 
present the result in two limiting cases. If the velocity is 
not too large, Vd << ceek^^Rp/l^ the main contribution 
to the k integral comes from the small k region, corre- 
sponding to length scales much larger than the range of 
the pinning potential. The upper limit of the wavevector 
integral can be extended to infinity, with the result, 



SvdiF) 



J2 nG)G'G ■ Fsgn(G • v^) 

47C66 
Tn,, 



(8) 



In this intermediate velocity regime the lifetime of elas- 
tic deformations of the sliding medium is small compared 
to the time to cross the range of the pinning potential, 
yielding collective pinning of the lattice. For the two- 
dimensional case considered here in this regime one ob- 
tains a force-independent correction to the drift velocity. 
Conversely, if Ud >> ceekg^Rp/j, the time needed to 
cross the range of the pinning potential - and therefore 
to see uncorrelated disorder - is short compared to the 
lifetime of elastic deformations connecting neighboring 
vortices, and one obtains single-particle response, with 
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(9) 



Notice that the coefficient of l/F in Eq. || is indeed 
independent of vortex density. 



IV. NUMERICAL RESULTS 

We have performed molecular dynamics simulations 
of arrays of 300, 920 and 1200 vortices using periodic 
boundary conditions. The results presented below are 
for two sets of parameter values, unless otherwise spec- 
ified. The data for the array with a dilute concentra- 
tion of pins {Np/Ny = 0.5) are obtained with Ny — 920, 
Np = 460 and Ry = 9.9. For these parameter values the 
clean Abrikosov lattice has lattice constant ao — 3.54, 
with UyR^ = 9.0, and shear modulus cqq = 0.271. As dis- 
cussed earlier, all lengths are measured in units of Rp and 
forces are in units of /„ . The data for the densely pinned 
array {Np/Ny — 133) are obtained with Ny — 300, 
Np = 40,000 and Ry = 20. In this case Oq = 7.44, 
with UyRl — 8.3, and cqq — 0.334. Table 1 shows the 
values of the collective pinning length Rc given in Eq. ^ 
and the threshold force estimated using the dimensional 
analysis discussed in section III for various pinning forces 



fp. For each set of parameters the value of the threshold 
force given in the table is the smaller of the two estimates 



and 



Table 1 



Np/Ny = 0.5 


Np/Ny = 133 


fp Rc/ Qq 


fp Rc/ OO Pj''^ 


0.2 20.0 2.1 X lO-'^ 


0.03 7.8 1.4 x 10--* 


0.5 13.3 1.3x10"^ 


0.1 2.3 1.6 X 10"'' 


1.0 4.0 5.3 X 10-^ 


0.3 0.8 1.0x10-^ 


1.5 2.7 1.2x10-^ 


3.0 0.08 1.0 x 10-^ 


2.0 2.0 2.1 X 10-^ 





The drift velocity of the vortex array is shown as a 
function of driving force in Figs. 2a and 3a for var- 
ious values of the maximum pinning force fp. Figure 
2 is for the sample with a low concentration of pins 
[Np/Ny ~ 0.5), while Fig. 3 is for the densely pinned 
sample {Np/Ny — 133). Figures 2b and 3b display the 
differential resistivity dvd/dF. Both Figs. 2a and 3a 
show a systematic evolution of the shape of the VI curve 
with increasing disorder strength not unlike that ob- 
served in the V-I curves of real superconductors Q| . For 
small pinning forces the velocity is nonlinear in F only 
very near threshold, where it exhibits a very small region 
of negative curvature. Correspondingly, there is no peak 
in the differential resistivity. At larger pinning forces 
there is a change in the sign of the curvature of the mean 
velocity that occurs at a value Fp^ak above threshold, but 
well in the nonlinear region, and yields a peak in the dif- 
ferential resistivity. The location of this peak moves to 
larger driving forces as the pinning force increases. This 
dependence is particularly strong in the sample with a 
dense pin array. For F > Fp^ak the VI characteristic is 
concave down and as F grows it approaches the asymp- 
totic value Vd = F, corresponding to a freely sliding ar- 
ray. In this region the deviations from the linear behav- 
ior Vd — F are fit quantitatively by the single-particle 
perturbation theory result given in Eq. |. This can 
be rewritten as Vd/F « 1— < Fp{0) > /F^, where 
< Fp{0) >— T/Rp is the mean square total pinning 
force. For the specific pinning potential used in our sim- 
ulations, < F^{0) (7r/30)npi?^j2^ xhe rms velocity 

fluctuations defined as Vrms =< [j{ X^i • F — Vd]'^ >^^^ 
are also fit quantitatively by their single particle value, 
Vrms = [< Fp > /2Ny]^^^ in this region. We stress that 
for very strong pinning the flux array can be very disor- 
dered even in the region F > Fp^ak, with sometimes as 
much as 50% of the vortices with a coordination number 
different from 6 (see Fig. 4b below). This is because dis- 
locations can be frozen in the sliding lattice in our T = 
simulations, yielding a disordered array that slides as a 
whole, with dislocations moving along at the same ve- 
locity as the rest of the lattice. This behavior may be 
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a finite-size effcet and is related to the hysteresis in the 
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FIG. 3. Drift velocity (a) and differential resistivity dva/dF 
(b) versus driving force for the array with a dense distribution 
of pinning centers, Np/Nv = 133. The curves obtained by 
FIG. 2. Drift velocity (a) and differential resistivity dv^/dF ramping the force up and down are indistinguishable. The 
(b) versus driving force for the array with a dilute distribution ^''''O'" ^ars represent the value of v^ms- Notice that both va 
of pinning centers, N^/N, = 0.5. The curves obtained by '^"'^ ^ ^ave been divided by ,U to display the data obtained 
ramping the force up and down are indistinguishable. The '^^^'^^^^ values of the pinning force on the same graph, 

error bars represent the value of Vrms- 

For Np/Ny = 0.5 the threshold force remains very 
small for all values of the pinning force studied. We 
cannot exclude that the threshold force may vanish in 
this case, but we have not performed extensive runs near 
threshold and finite-size scaling to determine the loca- 
tion of the threshold precisely. In fact, for a model with 
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dislocations, it is not necessary that Ft ^ at zero tem- 
perature (T = 0). In contrast, the threshold force is large 
and clearly nonzero for Np/Ny = 133. While we have not 
determined the threshold value accurately, we find that 
the numerical estimate agrees in order of magnitude with 
the dimensional estimates given in Table 1. Our numeri- 
cal V-I characteristics resemble those obtained in earlier 
simulations of the same model [[7|,p|JTl|] . There is, however, 
an important difference between our results and those of 
Koshelev and Vinokur in that we see no hysteresis in 
the V-I characteristics other than that arising from finite 
size effects. For small systems, we do see hysteresis in 
the V-I curves similar to that reported by Koshelev and 
Vinokur, but this hysteresis vanishes in larger systems. 
This hysteresis may be due to the periodic boundary con- 
ditions, which lead to metastable periodic attractors for 
the dynamics: the hysteresis in small systems in our sim- 
ulations is associated with such periodic attractors. 

We have also studied the evolution of the number of de- 
fects in the lattice with driving force. In most of our runs 
the flux array is prepared in an initial random disordered 
configuration and the driving force is then ramped up 
from zero. We have tracked the defects in the lattice by 
doing Voronoi constructions during the run and counting 
the number of vortices that are not 6-fold coordinated. 
The 5- and 7-fold coordinated vortices are disclinations 
in the two-dimensional triangular flux lattice and, when 
paired, correspond to a dislocation. Fig. 4(a) shows the 
time-averaged fraction of 6- fold coordinated vortices as a 
function of driving force for Np/N^ =0.5. In looking at 
these curves it should be kept in mind that the number 
of defects present at F = depends on initial conditions 
and the realization of disorder. As observed earlier [ pT| , 
we find that the flux array orders at large driving force. 
The value of driving force where the number of defects 
starts dropping is of the order of the location of the peak 
in the differential resistivity. Again, we have observed no 
hysteresis in the number of defects for the sample with 
a low concentration of pins other than one arising from 
finite size effects. We do, however, find hysteresis in the 
number of defects when the disorder is very strong, as 
shown in Fig. 4b. In this case when we ramp-up the 
force from an initial disordered vortex configuration, de- 
fects get "frozen in" and the lattice never orders. Starting 
from an ordered configuration at high fields, the lattice 
maintains its order. We have not been able to exclude 
that this hysteresis is also a finite-size effect. We ex- 
pect that this hysteresis will disappear in the presence of 
thermal fluctuations. The value of driving force where 
the defects become frozen-in apparently corresponds to 
the onset of the single-particle behavior discussed below. 
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FIG. 4. Fraction of 6-fold coordinated vortices versus driv- 
ing force. Fig. (a) is for Np/Nv — 0.5 and the parameter 
values of Fig. 2. In this case the curves obtained by ramping 
the force up and down are virtually indistinguishable and no 
hysteresis is observed. Fig. (b) is for Np/N^ = 133, with the 
parameter values of Fig. 3 and fp = 3. The lower curves 
are obtained by ramping up the force from an initial disor- 
dered configuration of the flux array. Data for both N^j = 300 
(circles) and = 1200 (triangles) are shown to display the 
finite size effect. The upper curve (square) is obtained by 
ramping down the force from an initial ordered configuration 
with = 300. 
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In order to correlate the macroscopic response of the 
vortex array with the details of the microscopic vortex 
motion, we have performed a variety of visualizations 
and we have studied the spatial distribution of vortex 
velocities. One method of displaying the data that we 
have found useful is to plot histograms of the component 
of the vortex velocity in the direction of the driving force 
(x-component) . Figure 5 shows the evolution of such his- 
tograms of the x-component of the instantaneous vortex 
velocity near threshold with driving force for parameter 
values that yield crinkle flow. For all driving forces, the 
histograms have a single maximum at a value of velocity 
close to the mean drift velocity. The location of the max- 
imum moves to larger velocities as the driving force in- 
creases. Very close to threshold velocities are distributed 
asymmetrically about the mean value and the histograms 
are not unlike those obtained from a phase-only model of 
CDW's As the driving force increases the velocity 

distribution becomes sharper and symmetric. 
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FIG. 5. Evolution of histograms of instantaneous velocity 
Vx with applied driving force for the parameters of Fig. 2 and 
fp = 0.2. The system exhibits crinkle flow near threshold. 

The instantaneous velocity distributions are quite dif- 
ferent for parameter values where plastic flow is found. 
In this case the velocity histograms display a clear bi- 
modal structure near threshold, as shown in Fig. 6, in- 
dicating the presence of two distinct "typical" velocities 
of the vortices. The first peak, which is approximately 
centered at zero very near threshold, is determined by 
the "slow" vortices in the array, located in pinned or 
temporarily pinned regions. This peak has finite width 
due to the oscillations of the "slow" or pinned vortices 
about the pins due to interactions with vortices flowing 



nearby. The peak at larger velocity is determined by the 
"fast" vortices that flow in channels around the pinned 
regions. We stress, however, that individual vortices are 
sometimes "slow", sometimes "fast". 
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FIG. 6. Evolution of the histograms of instantaneous ve- 
locity Vx with applied force for the parameter values of Fig. 
2 and fp — 1. The system exhibits plastic flow near threshold 
and the histograms clearly display a bimodal structure in this 
region. 

As a result of the temporal fluctuations in the ve- 
locity, we expect a large voltage noise in this region, 
perhaps related to the experiments by Marley, Bhat- 
tacharya and Higgins |^ . The various curves correspond 
to different driving forces, ranging from close to thresh- 
old (the threshold force for this case is estimated as 
Fx ~ 5.3x10""^) to well within the linear regime. The 
bimodal structure disappears at a value of F close to the 
location Fp^ak of the peak of the differential resistivity 
(here Fpeak ^ 0.125). Beyond this value the V-I curve 
rapidly becomes linear and the velocity distribution is 
narrow and symmetric, centered at the mean velocity. 
The origin of the maximum in the differential resistivity 
can easily be traced back to the shape of the velocity his- 
tograms by studying the location of the two peaks and 
their relative weights (Fig. 7) as functions of the driving 
force. Using a crude approximation, we can write the 
drift velocity of the vortex array as = UgVs + ^fVf, 
where and Uf are the fraction of slow and fast vor- 
tices, respectively, identified with the area under the re- 
spective peak of the velocity distribution and shown in 
Fig. 7(b). Similarly, and Vf are the corresponding ve- 
locities, identified here with the location of the two peaks 
and displayed in Fig. 8a. Using rt^ + nf = 1, we obtain 
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(dvd/dF) « {duf /dF){vf — Vg). For F < Fpeak the slow 
vortices are essentially not moving {vg ~ 0) while the 



0.40 r- 

0.35 - 



0.30 - 




0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 
F 



1 .0 




0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 
F 



FIG. 7. These figures display the evolution with applied 
force of the location (a) and weights (b) of the peaks of the 
histograms of Fig. 6. The location is obtained by recording 
the velocities at which each peak is maximum. The weight 
of the F** peak is the ratio of the area under the part of the 
histogram up to the minimum between the two peaks, to the 
total area under the curve. The weight of the 2"'' peak is 1 — 
the weight of the F* peak. 

number of fast vortices is increasing rapidly, leading 
to a superlinear V-L Above Fpeak the rate duf/dF at 



which the slow vortex fraction decreases and the slow 
vortex fraction increases slows down considerably, while 
(vf — Vs) ^ - it is this slowing down of the rate dn / /dF 
with increasing F that is responsible for the peak in the 
differential resistivity. 

The bimodal structure of the velocity histograms dis- 
cussed above reflects the spatial inhomogeneity of the 
instantaneous vortex velocity. A crucial question for the 
characterization of plastic flow is whether or not this bi- 
modal structure will persist when the vortex velocity is 
time averaged over times larger than those corresponding 
to an average displacement of the lattice of at least a lat- 
tice constant. It has been suggested that an important 
distinction between plastic and elastic response can be 
found in the correlations of the time-averaged velocity [|| . 
Indeed in a model where dislocations are forbidden and 
the response is therefore elastic, the time-averaged ve- 
locity will be spatially homogeneous and correlated over 
the entire system size. In contrast, in a system exhibiting 
plastic flow the time-averaged local velocity should still 
be spatial inhomogeneous, yielding bimodal structure of 
the corresponding histogram. We have constructed his- 
tograms of time-averaged vortex velocities, defined as 
Vi(T) — ^Vi{t) for various values of T, where T ~ 1 
yields the histogram of instantaneous velocity discussed 
above. The histograms are shown in Fig. 8 for two val- 
ues of the driving force. The bimodal structure clearly 
remains for the time-averaged velocities. 

To summarize, our model of a two-dimensional flux ar- 
ray driven through quenched disorder exhibits two types 
of response. For very weak disorder strength, the flux ar- 
ray exhibits "crinkle" motion, with correlated patches of 
vortices making small forward jumps at different times, 
like a tablecloth being pulled on a rough surface. At 
very small driving forces the lattice contains an appre- 
ciable number of 5- and 7- fold disclinations. The defects 
are concentrated at the boundaries between the patches 
with correlated velocities and their number drops rapidly 
to zero with increasing driving force. The distribution of 
vortex velocities displays a single maximum at a value of 
the velocity of the order of the mean velocity of the array 
(Fig. 3). As the driving force is increased and the flux 
array gets depinned, the maximum shifts to a higher ve- 
locity and the distribution broadens with no substantial 
change in shape. The histograms of time-averaged veloc- 
ities become sharper as the averaging time increases and 
stop changing once the averaging time exceeds the time 
over which the vortex lattice advances on the average a 
distance of a few Rp. This type of response is similar to 
that observed by Hu and Westervelt in magnetic bubble 
arrays These authors report observing a bimodal 

velocity distributions, but this is because they only look 
at the distribution of velocity over a very small time 
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FIG. 8. These figures compare histograms of instantaneous 
and time-averaged vortex velocities. The vortex velocities are 
averaged over a time T. The value T = 1 yields the in- 
stantaneous velocity histograms of Figs. 5 and 6. Figure 
(a) corresponds to parameter values yielding "crinkle" flow 
(same as Fig. 3, for fp = 0.03 and F = 0.0054). In this 
case the instantaneous velocity histogram exhibits a single 
maximum near Vd ~ .001 and is essentially identical to the 
histogram obtained with T = 10 and T = 160. Notice that in 
a time T = 10 the vortices displace on the average a distance 
Ax ~ O.OI800. The histogram for T = 640 is sharper, but 
qualitatively unchanged, indicating that the system motion 
is well correlated for these parameter values. The histograms 
shown in Fig. (b) refer to parameter values yielding plastic 
flow (same as Fig. 3 for fp = 3 and F = 3.5). In this case 
the histogram of instantaneous velocities (T = 1, solid curve) 
has a bimodal structure, which persists upon time averaging 
(T = 5 or Ax- ^ lAao and T = 10 or Ax ~ 2.7oo), indicating 
true plastic response. 

scale, smaller than the time required for the array as a 
whole to advance a distance of the order of the range 
of the pinning potential. Hu and Westervelt also argue 
that the response exhibits scaling in this regime. While 
we have not studied in sufficient detail the region near 
threshold to observe scaling, it seems quite plausible that 
this "crinkle" regime will exhibit generic critical behav- 
ior analogue to that predicted and observed for an elastic 
medium, in spite of the presence of defects. The fairly 
large number of defects present in our system at very 
small driving forces may very well be an artifact of our 
initial conditions, with vortex positions chosen at ran- 
dom, rather than equilibrated. It may be that if a low 
but finite temperature is introduced in the model and the 
flux array is initially allowed to equilibrate in the pres- 
ence of the disorder, the number of defects present for 
the parameter values yielding crinkle motion would be 
practically negligible even at the smallest driving forces. 

In most of the region of parameter space studied we 
have observed plastic flow of the flux array. This regime 
is characterized macroscopically by a change in the sign of 
the curvature of the V-I characteristic well above thresh- 
old, which yields a maximum in the differential resistance 
dvd/dF, and by a large number of defects in the region 
below the maximum. The flow is spatially inhomoge- 
neous. Over a short time interval one observes fluid-like 
flow of moving regions around pinned regions. On the av- 
erage, however, all vortices participate in the motion in 
the sliding state and no regions of the array are stuck for 
the entire length of the simulation, even near threshold. 
The evolution of the velocity distribution with driving 
force is shown in Fig. 4. There is clearly a bimodal ve- 
locity distribution in the sliding state which persists until 
all defects have healed and the V-I has become linear. A 
single vortex is in general "slow" for some of the time 
and "fast" for some of the time. This should result in a 
large voltage noise, as observed by Marley, Bhattacharya 
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and Higgins 

As the pinning force is increased, the persistence time 
of this structure of pinned and flowing regions grows 
and pinned vortices remain pinned for longer and longer 
times. For the situation of strongest pinning among those 
studied (/p = 3 for Np/Ny — 300), we find that some 
vortices are pinned for the entire length of the simula- 
tions, while others are flowing quite freely in channels 
surrounding the pinned regions. In this case the struc- 
ture of the channels is time independent near threshold. 
As the driving force is increased all vortices are eventu- 
ally depinned and the V-I becomes linear. The array is 
very defective near threshold and the velocity distribu- 
tion is bimodal and remains bimodal when velocities are 
time averaged. The flux array displays a filamentary mo- 
tion similar to that observed by Middleton and Wingreen 
for an array of quantum dots ||29[| , where near threshold 
the current flows in a single narrow channel and exhibits 
critical scaling. Critical scaling was also predicted in a 
model of fluid flow down a rough incline by Narayan and 
Fisher ^ . In this case again the flow pattern consists of 
directed channels that run across the system. It may be 
that in this very strong pinning regime the driven flux 
array also exhibits generic critical behavior, not unlike 
that of the fluid model of Ref. ^ . 

Some insight on the parameter regions where the two 
regimes described above may be expected to occur can 
be gained by considering the Larkin-Ovchinnikov pin- 
ning length, given in Table 1 for the case where the 
range Rp of the pinning potential is small compared 
to both Ry and oq. If Rv » clq, the shear modulus 
of the flux array can be estimated as Cgg ~ fvRv /^o 
and Lc ^ RvfvV^v/ {fp\/Np). The flux array should 
therefore be pinned collectively (Lc >> ao) provided 
fv V^v /{fp y/Np ) > 1. We find that these inequalities 
are generally satisfied in the cases where we observe crin- 
kle response. Conversely, when fvVN^/{fp\/^) << 1, 
even if i?^ >> oq, we obtain Lc ~ ao and vortices can be 
pinned individually yielding plastic flow or eventually fil- 
amentary flow. A transition between these two regimes is 
indeed observed as fp is increased for the parameter val- 
ues discussed above (here Lc decreases over almost two 
orders of magnitudes over the range of pinning forces 
studied, from Lc ^ 20ao for fp = 0.2 to Lc ^ O.Sao for 
fp = 10.). Using again the above estimate for cgg for 
Ry >> ao we find that in the collective pinning regime 
the threshold force needed to depin the array is given 
by Ft ~ {Np/Ny){Rp/Ry){fl/fy). The increase of the 
threshold force with fp observed in Fig. 2a is consistent 
with this dependence. 
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